function J_i = HJB_t( J_old, Q_t, param,aux)

    ind_stay = 1   + (aux.n_m+1).*(0:aux.n_m-1);

    %bounds for loop over m_h
    a= log(param.m_h)-0.05;
    b= log(param.m_h)+0.05;

    diff=1;
    while diff>aux.tol_HJB

        param.ln_m_HJB_r = (a+b)./2; % guess for m_h

        A= T_lam_matrix_HJB_t(param.ln_m_HJB_r,param.ln_m_HJB,Q_t,param,aux); %
        param.T_lam=A.T_lam;
        param.T_lam2=A.T_lam2;

        T_tr=param.T+param.T_lam.*(param.ln_m_HJB<param.ln_m_HJB_r)+param.T_lam2.*(param.ln_m_HJB<param.ln_m_HJB_r);
        A=1./aux.dt*speye(aux.n_m,aux.n_m)+aux.CN.*T_tr;
        B=1/aux.dt.*speye(aux.n_m,aux.n_m)-(1-aux.CN).*T_tr;

        J_i=J_old;
        J_i_old=J_i;

        diff_J=1;
        u=1;
        while diff_J>=10^(-10)

            income   = (1-param.omega_1).*exp(param.ln_m_HJB) -param.omega_0 - param.T_lam2*((param.ln_m_HJB>=param.ln_m_HJB_r).*min(param.c,J_old)) + B*J_old;

            P_m=aux.pen.*(J_i<0);
            P_p=aux.pen.*(J_i>param.c+param.K);%.*(param.c-J_i_old);

            AA=A;
            AA(ind_stay)=AA(ind_stay)'+P_m+P_p;
            J_i=AA\(income+P_p.*(param.c+param.K)+P_m.*(0));

            diff_J=max(abs(J_i_old-J_i));
            J_i_old=J_i;
            u=1+u;
            if u>10^4
               error('HJB failed to coverge')
            end
        end
        m_tmp=interp1(J_i_old,param.ln_m_HJB,param.c);
        diff=(m_tmp-param.ln_m_HJB_r)^2;
        if param.ln_m_HJB_r>m_tmp
            b=param.ln_m_HJB_r;
        else
            a=param.ln_m_HJB_r;
        end
    end
end

